C.................... MSIS00-SUBS.FOR ................................
C..... this program evaluates the solar zenith angle, determines the
C..... the day and time and calls Hedins msis model(s) to get
C..... the neutral densities and temperature.
      SUBROUTINE AMBS(J,DALT,GL,DD,T,SZA,SAT,GLON,GLATD)
      IMPLICIT REAL(A-H,O-Z)
      REAL NDFNOR,NDFSOU,TINFN,TINFS,DELTIN,DELTIS,TMSN,TMSS
      REAL EXFACN,EXFACS,DEC,DIP,HDFACN,HDFACS
      REAL RALT,GLATD,AE_OP,AE_OX,AE_N2,AE_O2,AE_NO,AE_TE,AE_TI
      DOUBLE PRECISION DALT,GLN,GL
      INTEGER TPRIN,TMODON
      !.. DD is used to transfer MSIS densities D out of AMBS so that there is
      !.. no need to expand the arrays outside those routines that call MSIS
      REAL T(2),D(9),DD(19)
      REAL D350M(9),D350S(9)  !.. used in modified MSIS
      !... TN,OX,O2,N2,HE,H perturbations from GWs, Nth, East, Ver, Bmer winds
      REAL GWDENTN(9),GWIND(3),GW_BMER    !.. Added 2023-02-06
      COMMON/MSIS/AP(7),DEC,ETRAN,BLON,F107,F107A,IDAY,SEC,GLN(401)
      COMMON/SOL/UVFAC(59),EUV
      COMMON/GPAR/RE,PLAT,PLON,PI,ANGVEL
      COMMON/NDFAC/TPRIN,NDFNOR,NDFSOU,TINFN,TINFS,DELTIN,DELTIS
     >   ,TMSN,TMSS,EXFACN,EXFACS,TMODON,HDFACN,HDFACS
      !.... common suggested by Kent Miller. zeroed in common block
      COMMON/AMBCOM/SECSAV,NY,NDY,LY,SX,ID,DECL
      COMMON/WPAR/NWTIMN,NWTIMS,NCOLN(9),NCOLS(9),UTOLTN,UTOLTS
     >  ,HMOWNN,HMOWNS

      !..... Since IDAY does not change and UT continually increases, need to
      !..... evaluate actual day (ID) and UT (SX) for MSIS.
      CALL ACTUAL_DAY(IDAY,SEC,ID,SX) 

      CALL SOLDEC(ID,SX/3600.,DECL,ETRAN)  !.. solar declination, ephemeris transit
      
      !....... transform magnetic to geographic coordinates---
      CALL BTOG(IEXP,BLON,SNGL(GL*57.296),GLATD,GLON)
      CALL BTOGAX(J,GLATD,GLON,DEC,DIP)   !.. Apex values
      GLAT=GLATD/57.29578

      !...... Evaluate Solar apparent time
      SAT=(SX-ETRAN+43200.0)/3600+GLON/15.0
      IF(SAT.LT.0) SAT=SAT+24
      IF(SAT.GT.24) SAT=SAT-24

      !...... Evaluate Solar zenith angle
      HH=(SAT-12.)*15.0*PI/180.       !.. hour angle
      !.. Evaluate argument for ACOS and test because roundoff error 
      !.. may cause invalid ACOS argument      
      SZA_ARG=COS(GLAT)*COS(DECL)*COS(HH)+SIN(GLAT)*SIN(DECL)
      IF(SZA_ARG.GT.0.999) SZA_ARG=0.999
      IF(SZA_ARG.LT.-0.999) SZA_ARG=-0.999
      SZA=ACOS(SZA_ARG)                 !.. solar zenith angle

c...... For studying the effects of nightime EUV
c...      IF(ABS(GLATD-38).LT.3.AND.SZA.GT.93.0/57.29578) SZA=93.0/57.29578
c...      IF(ABS(GLATD+35).LT.3.AND.SZA.GT.95.0/57.29578) SZA=95.0/57.29578
      
      RALT=SNGL(DALT) !.. convert double to single precision for MSIS

      !.. This section is used with FLIPTN to keep the MSIS O density. See also
      !.. DOX statement further down
      !.. DATA DOX/-9/   !.. initialize normal O density
      !.. D(1)=-2
      !.. This call to get normal MSIS densities and temperatures
      !.. CALL GTD7(ID,SX,RALT,GLATD,GLON,SAT,F107A,F107,AP,48,D,T)
      !.. DOX=D(2)    !.. normal O density
                
      !.. TMODON is set in ALTWIN and passed to MSIS to switch on greater Kp response
      D(1)=TMODON
      !.. This switches on the PGR modification to MSIS not based on NmF2 if = -1
      IF(UVFAC(59).EQ.-2) D(1)=0   !-1

      !.. This call to get MSIS densities and temperatures
      CALL GTD7(ID,SX,RALT,GLATD,GLON,SAT,F107A,F107,AP,48,D,T)

      !.. The following routine allows the change of individual MSIS parameters
      !..CALL MODMSIS(ID,SX,RALT,GLATD,GLON,SAT,F107A,F107,AP,D,T)

      !.. The following routine allows simulation of hot O
      !..CALL OXHOT(ID,SX,RALT,GLATD,GLON,SAT,AP,F107A,F107,D,T) !.. NOT tested
      
      !.. Transfer densities D to DD which can be scaled through UVFAC array
      !.. Species order = 1-He, 2-O, 3-N2, 4-O2, 5-Ar, 7-H, 8-N
      !.. D(6) = total mass density. D(9) = Hot O density
      DO 50 I=1,5
         DD(I)=D(I)*ABS(UVFAC(50+I))
 50   CONTINUE
      !.. FLIPRUN.DDD H density multiply factors for both hemispheres 2022-05-29
      DD(7)=D(7)*ABS(UVFAC(50+6))  !.. Use same factor in both hemispheres
      !.. Next line gives changes the H factor in the south if desired
      IF(GLATD.LT.0.0.AND.UVFAC(46).GT.0.01) DD(7)=D(7)*ABS(UVFAC(46))    !.. Sth

      !.. H multiply factors from column 7 of the hmF2/wind data file 2022-08-15
      IF(GLATD.GE.0.0.AND.ABS(HDFACN-1.0).GT.0.001) DD(7)=D(7)*HDFACN     !.. Nth
      IF(GLATD.LT.0.0.AND.ABS(HDFACS-1.0).GT.0.001) DD(7)=D(7)*HDFACS     !.. Sth
 
      DD(8)=D(8)*ABS(UVFAC(50+7))  !... N density
      DD(9)=D(9)                   !... Hot O density
      
      !... Get Gravity Wave data to modify neutral densities and wind
      !... TN,OX,O2,N2,HE,H perturbations from GWs. GW added 2023-02-06
      DO I=1,9
         GWDENTN(I)=0.0
      ENDDO
      CALL GRAVITY_WAVE_DATA
     >    (0,J,SEC/3600.0,RALT,GLATD,GWDENTN,GWIND,GW_BMER)
         DD(1) = DD(1)+ GWDENTN(5)   !...  He
         DD(2) = DD(2)+ GWDENTN(2)   !...  OX
         DD(3) = DD(3)+ GWDENTN(4)   !...  N2
         DD(4) = DD(4)+ GWDENTN(3)   !...  O2
         DD(7) = DD(7)+ GWDENTN(6)   !...  H
         T(2) =  T(2) + GWDENTN(1)   !...  Tn

      !-- MSIS is being adjusted to improve NmF2
      !-- scale [O] and [N2] to improve model NmF2. See ALTWIN
      IF(GLATD.GE.0) THEN
          IF(NDFNOR.LT.0.3.OR.NDFNOR.GT.2.0) WRITE(6,*) ' ** in AMBS'
          DD(2)=DD(2)*NDFNOR**1.000
          DD(3)=DD(3)/NDFNOR            !.. Change N2
          DD(4)=DD(4)/NDFNOR            !.. Change O2
          !..DD(7)=DD(7)*NDFNOR         !.. Change [H]
      ELSE
          IF(NDFSOU.LT.0.3.OR.NDFSOU.GT.2.0) WRITE(6,*) ' ** in AMBS'
          DD(2)=DD(2)*NDFSOU**1.000     !.. Change O
          DD(3)=DD(3)/NDFSOU
          DD(4)=DD(4)/NDFSOU
          !..DD(7)=DD(7)*NDFSOU
      ENDIF

      !.. IF(DOX.GT.0.0) DD(2)=DOX  !.. Set O density to normal MSIS

      !.. This section allows FLIP to substitute measured values
      !.. from a file. If the data are bad,
      !.. zero values are returned and FLIP values are retained.
      CALL GET_AE_DENS(RALT,SEC/3600.,GLATD,AE_OP,AE_OX,AE_N2
     >    ,AE_O2,AE_NO,AE_TE,AE_TI)
      IF(AE_OX.GT.100) DD(2)=AE_OX
      IF(AE_N2.GT.100) DD(3)=AE_N2
      IF(AE_O2.GT.100) DD(4)=AE_O2

      RETURN
      END
C:::::::::::::::::::::::::::: MODMSIS :::::::::::::::::::::::::::::::::::::
C... This routine is used to modify individual MSIS densities and temperatures
C... to study the effects 
C... Written by P. Richards in December 2002
      SUBROUTINE MODMSIS(ID,SX,RALT,GLATD,GLON,SAT,F107A,F107,AP,D,T)
      IMPLICIT NONE
      INTEGER ID,ISW,I
      REAL SX,RALT,GLATD,GLON,SAT,F107A,F107,AP(7),D(9),DP(9),T(2),TT(2)
      REAL APP(7),SW,SWC,JTI
      COMMON/CSW/SW(25),ISW,SWC(25)
      DATA JTI/0/
      JTI=JTI+1
      IF(JTI.EQ.1) WRITE(6,*) ' ~~~~ MODMSIS is CALLED ~~~~~'

      !.. Set baseline thermosphere using constant AP      
        DATA APP/7*15/
      SW(9)= 1.0        !.... allow MSIS daily AP
      D(1)=1
      CALL GTD7(ID,SX,RALT,GLATD,GLON,SAT,F107A,F107,APP,48,D,T)

      IF(APP(1).GT.0) RETURN

      SW(9)= -1.0      !..  allow MSIS to use the 3 hour AP
      !.. This call to change one or more of the densities or temperatures
      !.. to the actual MSIS value. 

      !..DP(1)=-1    !.. switches on PGR enhanced MSIS
      CALL GTD7(ID,SX,RALT,GLATD,GLON,SAT,F107A,F107,AP,48,DP,TT)
       !D(2)=DP(2)     !.. O density
       !D(3)=DP(3)     !.. N2 density
       D(4)=DP(4)     !.. O2 density
       !T(1)=TT(1)     !.. Tinf
       !T(2)=TT(2)     !.. Tn
      RETURN
      END
C:::::::::::::::::::::::::::: OXHOT :::::::::::::::::::::::::::::::::::::
C... This routine is used to simulate hot O using the He scale height
C... Written by P. Richards in December 2002. (NOT tested)
      SUBROUTINE OXHOT(ID,SX,ALT,GLATD,GLON,SAT,F107A,F107,AP,D,T)
      IMPLICIT NONE
      INTEGER ID,ISW
      REAL SX,ALT,GLATD,GLON,SAT,F107A,F107,AP(7),D(9),DP(9),T(2),TT(2)
      REAL HOTO,HE_TO_O,DHOTO

      !.. Start hot O procedure
      DATA HOTO/0.0/ !.. ratio of hot O to cold O at 500 km
      !.. Call to MSIS to add in hot O. Ratio is determined at 500 km
      CALL GTD7(ID,SX,500.0,GLATD,GLON,SAT,F107A,F107,AP,48,DP,TT)
      HE_TO_O=DP(1)/DP(2)       !.. HE to O ratio at 500 km

      DHOTO=D(1)*HOTO/HE_TO_O
            
      IF(HOTO.GT.0.0001.AND.GLATD.GT.0.AND.ALT.GT.490.AND.ALT.LT.510.) 
     >   WRITE(6,541)  NINT(ALT),D(2),DHOTO,D(1),HE_TO_O
 541  FORMAT(' HOT O is on: ',I7,1P,9E10.2) 
      D(2)=D(2)+DHOTO      !.. add hot O

      RETURN
      END
